Difference between raster and vector representations

There are different ways to plot maps in R. The packages you use will depend on which your aim is and on the format of your data. Sometimes, geographical data can be simply in a .txt or .csv file, where you have one column for latitude and one for longitude and other variables or measurements corresponding to that location in the following columns. Yo can also have shape files, which are special files for geographical information. This files may contain points, lines or polygons. They are usually .shp files. Finally you can have raster files which are something similar to pictures or .png files because these files have pixels and (at least) a value associated with each pixel. The difference between a .png file and a shape file is that the shape-file has associated coordinates that locate this pixels on a specific surface of the world. The format is usually .geotiff or just .tiff.

In this meetup we will cover the basis of working with this formats. But this is only the starting point. There is a lot out there.

from link

Simplest example

We can just use geom_polygon from the library ggplot to create our first figure. LetÂ´s suppose we know the latitude and the longitude of all the nodes in our polygon. How can we plot it?

require(ggplot2)

funny <- data.frame(lat = c(62, 55, 48, 48, 62), long = c(20,10.5,20,10,10)) # from top right 
#following the clock
ggplot(funny, aes(x = long, y = lat)) + # we specify the data
  geom_polygon(fill="green") + # we plot it
  geom_point(aes(x=13.25, y=50.5), color="violet", size=18) # we can also add points this was

Information already availabe in R

The geometry of the polygons of the countries in the world are already loaded into R. You can access this info using the function map_data from the package ggplot2.

How do we plot only one or a subset of countries?

We need to get the data for the geometry of the countries. To do that we use the library maps and ggplot2.

countries <- c("Germany", "France")
# Get the data
some.maps <- map_data("world", region = countries)  # This function retrieves the data. The data we get is a df with lat and long of each node. 
head(some.maps)

##       long      lat group order  region subregion
## 1 14.21367 53.87075     1     1 Germany    Usedom
## 2 14.17217 53.87437     1     2 Germany    Usedom
## 3 14.04834 53.86309     1     3 Germany    Usedom
## 4 13.92578 53.87905     1     4 Germany    Usedom
## 5 13.90215 53.93896     1     5 Germany    Usedom
## 6 13.92168 53.99663     1     6 Germany    Usedom

require(dplyr)
require(tidyr)

# Mean of lat and long to then write the labels
lab.data <- some.maps %>% group_by(region) %>% summarise(long = mean(long), 
    lat = mean(lat))
lab.data

## # A tibble: 2 x 3
##   region   long   lat
##   <chr>   <dbl> <dbl>
## 1 France   3.23  46.2
## 2 Germany 10.4   51.2

Themes in ggplot In ggplot2 the display of non-data components is controlled by the theme system. The theme system of ggplot2 allows the manipulation of titles, labels, legends, grid lines and backgrounds. There are various build-in themes available that already have an all-around consistent style. A very common one for plotting data is theme_bw() and theme_map(). We use theme_map() which will give our plot the appearance of a map. To use theme_map() we need to load the package ggthemes().

require(ggthemes)

fra_and_ger <- ggplot(some.maps, aes(x = long, y = lat)) +
  geom_polygon(aes( group = group, fill = region))+ #plot polygon, use group to group all the nodes of the same country together. 
  geom_text(aes(label = region), data = lab.data,  size = 3, hjust = 0.5)

#using theme_map
fra_and_ger + #print labels
  theme_map() +
  theme(legend.position = "none")

#Using theme bw

fra_and_ger + #print labels
  theme_bw() +
  theme(legend.position = "none")+
  xlab("Longitude") + ylab("Latitude")

From a couple of countries to the whole world: Plotting simple and nice maps of the world

ggplot2 is becoming the standard for R graphs. However, it does not handle spatial data specifically. Handling spatial objects in R relies on Spatial classes defined in the package sp or sf. ggplot2 allows the use of spatial data through the package sf. sf elements can be added as layers in a graph. The combination of ggplot2 and sf enables to create maps, using the grammar of graphics,but incorporation geographical info.

1. We can just create a map from the world. We will do this using the library `ggplot2`.

The function ne_countries also retrieves info about the countries. But now you do not get just nodes with coordinates but an sf object that contains much more info.

require(rnaturalearth)
require(rnaturalearthdata)
require(sf)

world <- ne_countries(scale = "medium", returnclass = "sf")  # this is another function to get polygons of countries. 

class(world)

## [1] "sf"         "data.frame"

dim(world)

## [1] 241  64

head(world[, 1:5])

## Simple feature collection with 6 features and 5 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -70.06611 ymin: -18.01973 xmax: 74.89131 ymax: 60.40581
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##   scalerank      featurecla labelrank     sovereignt sov_a3
## 0         3 Admin-0 country         5    Netherlands    NL1
## 1         1 Admin-0 country         3    Afghanistan    AFG
## 2         1 Admin-0 country         3         Angola    AGO
## 3         1 Admin-0 country         6 United Kingdom    GB1
## 4         1 Admin-0 country         6        Albania    ALB
## 5         3 Admin-0 country         6        Finland    FI1
##                         geometry
## 0 MULTIPOLYGON (((-69.89912 1...
## 1 MULTIPOLYGON (((74.89131 37...
## 2 MULTIPOLYGON (((14.19082 -5...
## 3 MULTIPOLYGON (((-63.00122 1...
## 4 MULTIPOLYGON (((20.06396 42...
## 5 MULTIPOLYGON (((20.61133 60...

We can plot this geographical object using ggplot2 and geometry geom_sf

plot_w1 <- ggplot(data = world) + theme_bw() + geom_sf() + xlab("Longitude") + 
    ylab("Latitude") + ggtitle("World map", subtitle = paste0("(", dim(world)[1], 
    " countries)"))

plot_w1

This really looks like a map!

2. We can change the colors, add text …

plot_w1 +
   geom_sf(color="blue", fill="black") +
  theme(panel.background = element_rect(fill = "black"))+# Modify the theme to change the background 
  #Where was the funny polygon located?
  geom_polygon(data=funny, aes(x = long, y = lat), color="green", fill="green") + # we plot it
  # Add points using lat-lon information
  geom_point(aes(x=13.25, y=50.5), color="violet", size=2)

However our nice modern-art-style polygon looks different in the map compared to the first plot. Any idea about what is happening?

Coordinate reference systems (crs). Going from the rear world (the earth) to the simplifyed model (the map)

from link

After choosing elipsoid and datum we have to project from 3D to 2D.

from link

You can get a feeling of how the Mercator projection distorts our worldview at link.

Google and many apps use unprojected coordinates. When the coordinates are unprojected, they are in degrees and give the position on an sphere and not in a 2D surface. But, you still need an ellipsoid and datum to make clear which reference system you are using. All this information is coded in an EPSG code. The most common coordinate reference system crs (used by Google and most apps) is EPGS: 4326. When you have lat Lon coordinates and no more info, it is likely that the EPSG is 4326, which uses the datum WGS84.

To get a taste: EPSG: 4326 Coordinate Systems Worldwide

To find the one used in your region of interest:

Argentina Germany Europe

Make Choropleth Map

What is a Chorophlet Map? “A choropleth map is a thematic map in which areas are shaded or patterned in proportion to the measurement of the statistical variable being displayed on the map, such as population density or per-capita income.”(Wikipedia)

1. Load Data

require(readr)  # to use the function read_csv()
life.exp <- read_csv("data/LifeExp.csv")

The data was obtained from link, lot of cool data-sets are available for free.

2. Tidy the data set to have one column with country and another with life expectancy

life_exp <- life.exp %>%
  filter(year == 2016 & sex == "Both sexes") %>%  # Keep data for 2016 and for both sex
  dplyr::select(country, value) %>%                      # Select the two columns of interest
  rename(name = country, lifeExp = value) %>%   # Rename columns
  #We have a very common proble when using different sources, the name of the countries not allways match
  # Replace "United States of America" by USA in the region column
  mutate(
    name = ifelse(name == "United States of America", "United States", name), 
    name = ifelse(name == "Russian Federation", "Russia", name)
    )

3. Get data of the polygons of each country and merge map and life expectancy data

# attribute names
colnames(world)

##  [1] "scalerank"  "featurecla" "labelrank"  "sovereignt" "sov_a3"    
##  [6] "adm0_dif"   "level"      "type"       "admin"      "adm0_a3"   
## [11] "geou_dif"   "geounit"    "gu_a3"      "su_dif"     "subunit"   
## [16] "su_a3"      "brk_diff"   "name"       "name_long"  "brk_a3"    
## [21] "brk_name"   "brk_group"  "abbrev"     "postal"     "formal_en" 
## [26] "formal_fr"  "note_adm0"  "note_brk"   "name_sort"  "name_alt"  
## [31] "mapcolor7"  "mapcolor8"  "mapcolor9"  "mapcolor13" "pop_est"   
## [36] "gdp_md_est" "pop_year"   "lastcensus" "gdp_year"   "economy"   
## [41] "income_grp" "wikipedia"  "fips_10"    "iso_a2"     "iso_a3"    
## [46] "iso_n3"     "un_a3"      "wb_a2"      "wb_a3"      "woe_id"    
## [51] "adm0_a3_is" "adm0_a3_us" "adm0_a3_un" "adm0_a3_wb" "continent" 
## [56] "region_un"  "subregion"  "region_wb"  "name_len"   "long_len"  
## [61] "abbrev_len" "tiny"       "homepart"   "geometry"

life_exp2 <- left_join(world, life_exp, by = "name")

4. Plot with ggpoligon

life_exp_map <- ggplot(data = life_exp2) + geom_sf(aes(fill = lifeExp)) + scale_fill_viridis_c(option = "plasma", 
    trans = "sqrt")  # this allows you to choose different colour scale

life_exp_map

Clearly we do not have information for all the countries regarding life expectancy and/or we have some problems with non-matching country names. When we call left_join we only keep in the table the countries that are in the left table.

5. Add points

We can also add points to our plot. For example we could add some capitals. Points where samples were obtained, etc. We can also control the size and color of the points using another variable. For example, we could plot sample points and the size would be number of observations. We could plot one point per city and the size number of inhabitants…and so son and so forth.

country_capitals <- read_csv("data/country-capitals.csv")

south_america_capitals <- country_capitals %>% dplyr::filter(ContinentName == "South America")


south_america_capitals$CapitalLatitude <- as.numeric(south_america_capitals$CapitalLatitude)

require(ggrepel)

life_exp_map +
  geom_point(data=south_america_capitals, 
             aes(x=CapitalLongitude, y=CapitalLatitude),
             alpha=0.5)+
  geom_text_repel(data=south_america_capitals, # geom_text_repel avoidsoverlapping
            aes( x=CapitalLongitude, y=CapitalLatitude, label=CapitalName), 
            hjust=0, vjust=0, size= 3)+ 
  theme_bw()

6. Add scale and north

require(ggspatial)

life_exp_map + 
  annotation_scale(location = "bl", width_hint = 0.5) + # scale
    annotation_north_arrow(location = "bl", which_north = "true", #north arrow
        pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"),
        style = north_arrow_fancy_orienteering) +
    coord_sf(xlim = c(-20.15, 50.12), ylim = c(35.65, -45)) #crop

Rasters and polygons

Usually spatial data comes as a raster or as a polygon format. We usually have shape files with the onformation of sampling sites, neighborhoods, streets, rivers, hospitals, surveys, etc. saved as .shp files. Or we have lat-Lon data that we can transform into an sf object as we did before. We can also get raster files such as climate raster files or remote sensing data, land use, etc. Most of the times we want to plot all this info together (raster + one or more shape-files).

1. Load required libraries

require(rgdal)
require(raster)  #Package to work with raster files
require(sf)  #Package to work with shape files

All this libraries allow us to work with spatial data using R.

rgdal is a library that that allows us to read and write geospatial data in R. The library just translates the already existing library gdal into R. The link to gdalproject is link.
raster is the library that we use to work with raster formats.
sf encodes spatial vector data. We already used it.

2. Load the raster files and do some calculations

LetÂ´s suppose he have different sampling points in Germany. We are interested in the temperature difference between summer and winter in this sites. Can we do a map to display this info? How would you do that?

This are some usefull links to find data:

Link to the Federal Agency for Cartography Link to open BW data Link to WoldClim link to land use data Another useful link

All the data used in this example was download from this sites or directly from R.

t_july <- raster("data/wc2.0_5m_tavg_07.tif")  # loads the raster
t_december <- raster("data/wc2.0_5m_tavg_01.tif")  # loads the raster

# Get temperature range
t_diff <- abs(t_july - t_december)
plot(t_diff, main = "Temperature range")

Raster calculation is relatively straight forward using R if the raster have the same resolution. If the raster come from different sources you usually have to re-shape one raster.

more info here

3. Load the shape files

sites <- st_read("data/sample_sites.shp")  # loads the vector data - sites

## Reading layer `sample_sites' from data source `C:\Users\ichas\Documents\RLadies\geoR\rmarkdown\data\sample_sites.shp' using driver `ESRI Shapefile'
## Simple feature collection with 7 features and 3 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 8.081778 ymin: 47.73776 xmax: 13.30927 ymax: 53.69961
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

# get germany voundaries from world R data
germany <- world %>% dplyr::filter(name == "Germany")
plot(germany)

As you can see here, each polling has associated data. Geometry has the info to build the polygon. level, type, etc. is the info related to each polygon. If the Polygons have different info in this fields, the polygons will be displayed in different colors.

# get germany voundaries from world R data
gfs <- world %>% dplyr::filter(name %in% c("Germany", "Austria"))
plot(gfs)

4. Check reference system

We have data coming from different sources, we have to check that the coordinate reference system match (otherwise we will have lots of problems and we will be doing everything wrong)

# The function to get the coor. system is different between shapes and
# rasters
st_crs(germany)  #get the projection

## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"

st_crs(sites)  #get the projection

## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"

crs(t_diff)  #get the projection

## CRS arguments:
##  +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

st_crs(germany)$proj4string == st_crs(sites)$proj4string

## [1] TRUE

st_crs(germany)$proj4string == crs(t_diff, asText = T)

## [1] FALSE

The data do not have the same coord. reference system. We have to transform one.

sites_transformed <- st_transform(sites, crs = crs(t_diff, asText = T))  #transform coordinate sistem
germany_transformed <- st_transform(germany, crs = crs(t_diff, asText = T))  #transform coordinate

5. Crop the raster to the shape file extent

We get the extent (the coordinates of the extreme points) for our shape file of Germany.

germany_extent <- extent(germany_transformed)
germany_extent

## class       : Extent 
## xmin        : 5.85752 
## xmax        : 15.0166 
## ymin        : 47.27881 
## ymax        : 55.05874

# crop the land use data by the extend of BW. Crop funciton from the raster
# package
crop_tdiff <- crop(t_diff, germany_extent, snap = "out")
plot(crop_tdiff)

In order to get a raster that matches exactly with the shape, we have to mask it. We could also do this step without cropping first, but that would be much more computational intensive.

mask_tdiff <- mask(crop_tdiff, germany_transformed)
plot(mask_tdiff)

6. Plot with `tmap`

The library tmap is another grate option to build thematic maps. I find it very useful to work with raster data. Useful info can be found in the tmap vignette. Another library to plot raster data is rasterVis. We will not cover it today but you can check the vignette You can plot many layers together. Each needs some reference(point, raster, text, etc.)

require(tmap)

tmap_mode("plot") # static plot or leaflet. 

tm_shape(mask_tdiff) + 
  tm_raster (col= "layer" , style = "cont", n=10,  title = "T. Difference") +
tm_shape(germany_transformed) + #add ploygon of germany
  tm_borders() +
tm_shape(sites_transformed) +
   tm_symbols(col = "gray", scale = .5)+ ##add points
      tm_text("Site", size = 0.8,  ymod=-0.6, root=1, size.lowerbound = .60, #add text
        bg.color="gray", bg.alpha = .5) +
tm_layout(inner.margins = c(0.15, 0.30, 0.10, 0.1)) + # crop the extent of the map
tm_layout(legend.position = c("left","bottom")) + # add legent
    tm_compass() + # add north
  tm_scale_bar() #add scale

tm <- tm_shape(mask_tdiff) + 
  tm_raster (col= "layer" , style = "cont", n=10,  title = "T. Difference") +
tm_shape(germany_transformed) + #add ploygon of germany
  tm_borders() +
tm_shape(sites_transformed) +
   tm_symbols(col = "gray", scale = .5)+ ##add points
      tm_text("Site", size = 0.8,  ymod=-0.6, root=1, size.lowerbound = .60, #add text
        bg.color="gray", bg.alpha = .5) +
tm_layout(inner.margins = c(0.15, 0.30, 0.10, 0.1)) + # crop the extent of the map
tm_layout(legend.position = c("left","bottom")) + # add legent
    tm_compass() + # add north
  tm_scale_bar() #add scale


tmap_save(tm, filename = "t_range.png")

With tmap is really easy to do interactive maps!

require(tmap)

tmap_mode("view")

tm_shape(mask_tdiff) + tm_raster(col = "layer", style = "cont", n = 10, title = "T. Difference") + 
    tm_shape(germany_transformed) + tm_borders() + tm_shape(sites_transformed) + 
    tm_symbols(col = "gray", scale = 0.5) + tm_text("Site", size = 0.8, ymod = -0.6, 
    root = 1, size.lowerbound = 0.6, bg.color = "gray", bg.alpha = 0.5) + tm_layout(inner.margins = c(0.1, 
    0.3, 0.1, 0.1)) + tm_layout(legend.position = c("left", "bottom"))

## Reclasify values in the rasta flie. Reclasify function of raster file
## function

tmap is a grate library to make maps. Check the vignette

Hands on work!

Exercise 1: Geo-information available in R.

Get eh data for the following countries (or a subset) :Switzerland, Germany, Austria, Belgium, Netherlands, Denmark, Poland, Italy, Croatia, Slovenia, Hungary, Slovakia,Czech republic.
Plot each country with a different color.
Add the name in the center of the country.

Option 1 - Using only ggplot2

require(ggplot2)
require(tidyr)
require(dplyr)

# Some EU Contries
some_eu_countries <- c( "Switzerland", "Germany",
  "Austria", "Belgium", "Netherlands",
  "Denmark", "Poland", "Italy", 
  "Croatia", "Slovenia", "Hungary", "Slovakia",
  "Czech republic"
)
# Retrievethe map data
some_eu_maps <- map_data("world", region = some_eu_countries)

# Compute the centroid as the mean longitude and lattitude
# Used as label coordinate for country's names
region_lab_data <- some_eu_maps %>%
  group_by(region) %>%
  summarise(long = mean(long), lat = mean(lat))

ggplot(some_eu_maps, aes(x = long, y = lat)) +
  geom_polygon(aes( group = group, fill = region)) +
  geom_text(aes(label = region), data = region_lab_data,  size = 3, hjust = 0.5) +
  scale_fill_viridis_d() +
  theme_void()+ # This is what is defining the background and the colors
  theme(legend.position = "none")

Option 2 - Usinng sf objects

require(rnaturalearth)
require(rnaturalearthdata)
require(sf)
require(ggplot2)
require(ggthemes)
require(ggspatial)



world <- ne_countries(scale = "medium", returnclass = "sf")  # this is another function to get polygons of countries. 

# The names of the countries are different between the sourcecs. To check
# how to spell the names:

head(unique(world$name))

## [1] "Aruba"       "Afghanistan" "Angola"      "Anguilla"    "Albania"    
## [6] "Aland"

# Czech republic is averebiated in this case, we change it:

# Some EU Contries
some_eu_countries <- c("Switzerland", "Germany", "Austria", "Belgium", "Netherlands", 
    "Denmark", "Poland", "Italy", "Croatia", "Slovenia", "Hungary", "Slovakia", 
    "Czech Rep.")

some_eu_sf <- world %>% dplyr::filter(name %in% some_eu_countries)

# This is a possible solution
ggplot(data = some_eu_sf) + geom_sf(aes(fill = name)) + coord_sf(xlim = c(0, 
    25.12), ylim = c(35.65, 60)) + theme_bw()

# To get the names cetered

eu_points <- st_centroid(some_eu_sf)
eu_points <- cbind(eu_points, st_coordinates(st_centroid(eu_points$geometry)))

ggplot(data = some_eu_sf) + geom_sf(aes(fill = name)) + coord_sf(xlim = c(0, 
    25.12), ylim = c(35.65, 60)) + theme_map() + theme_bw() + geom_text(data = eu_points, 
    aes(x = X, y = Y, label = name), fontface = "bold", check_overlap = FALSE, 
    size = 3, color = "darkred") + theme(legend.position = "none") + xlab("Longitude") + 
    ylab("Latitude") + scale_fill_viridis_d() + annotation_scale(location = "bl", 
    width_hint = 0.3) + annotation_north_arrow(location = "bl", which_north = "true", 
    pad_x = unit(0.15, "in"), pad_y = unit(0.3, "in"), style = north_arrow_fancy_orienteering)

Option 3 tmap

require(tmap)

new_bb = c(-10, 35, 35, 60)
names(new_bb) = c("xmin", "ymin", "xmax", "ymax")
attr(new_bb, "class") = "bbox"

attr(st_geometry(some_eu_sf), "bbox") = new_bb


tmap_mode("plot")

tm_shape(some_eu_sf) + tm_borders() + tm_polygons("name") + tm_shape(eu_points) + 
    tm_symbols(col = "gray", scale = 0.5) + tm_text("name", size = 0.8, ymod = -0.6, 
    root = 1, size.lowerbound = 0.6, bg.alpha = 0.5) + tm_compass() + tm_scale_bar() + 
    tm_legend(show = FALSE)

Exercise 2: Geo-information from a CSV. Investigate the data frame. What can we learn from the data? which map could we do?

This data frame has information about all the R Ladies chapters around the world! It is quite interesting. There are many questions we can answer. But first we have to take a look at the data.

require(readr)


rladies <- read_csv("data/RLadiesChapters.csv")
head(rladies)

## # A tibble: 6 x 8
##      X1 chapter      date       members   lat     lon city         country
##   <int> <chr>        <chr>        <int> <dbl>   <dbl> <chr>        <chr>  
## 1     1 R-Ladies Ba~ 22/10/201~     389  41.4    2.17 Barcelona    ES     
## 2     2 R-Ladies Us~ 09/05/201~      23 -54.8  -68.3  Ushuaia      AR     
## 3     3 R-Ladies Bi~ 27/02/201~      38  43.2   -2.93 Bilbao       ES     
## 4     4 R-Ladies Ba~ 10/05/201~     114 -41.1  -71.3  San Carlos ~ AR     
## 5     5 R-Ladies Me~ 02/09/201~    1213 -37.8  145.   Melbourne    AU     
## 6     6 R-Ladies No~ 24/05/201~      46  45.2   19.8  Novi Sad     RS

# date stands for the time when the chapter was created. We will change the
# name of this column
rladies$created_at <- as.character(rladies$date)

Which would be a nice plot to show this data? Would you use a clorophlet map? Or Would you add points? There is no right or wrong! LetÂ´s try:

Option 1 – Only ggplot

Where are the chapters of R Ladies? Where are the bigger/smaller chapters?

require(ggplot2)
require(ggthemes)

# Here we do the same map we did at the begining of the Meetup
world <- ggplot() + borders("world", colour = "gray85", fill = "gray80") + theme_map()


map <- world + # We create points
geom_point(aes(x = lon, y = lat, size = members), data = rladies, colour = "purple", 
    alpha = 0.5) + scale_size_continuous(range = c(1, 8), breaks = c(250, 500, 
    750, 1000)) + labs(size = "Followers") + ggtitle("R-Ladies chapters around the world")

map

In which countries are there more R Ladies chapters? A clorophlet map.

require(dplyr)
require(tidyr)
# We summarise the number of chapters per country
country.chapters <- rladies %>% group_by(country) %>% summarise(n.chapters = n())

# We have to add this data to the world_map data that has the poligons

world_map <- map_data("world")

However we have a problem here…we have to join the data-sets but the country names are codded in a different way. This kind of problems are really common when we work with countries. I just found a txt with country name and code online here

countries <- read.csv("data/countries.txt")
colnames(countries)[2] <- "country"  #use the same name as in rladies df

rladies.country <- left_join(country.chapters, countries, by = "country")  # join rladies df with the names of the countries
colnames(rladies.country)[3] <- "region"

map.rladies <- full_join(rladies.country, world_map, by = "region")  # join the new rlades df with the map df

ggplot(map.rladies, aes(long, lat, group = group)) + geom_polygon(aes(fill = n.chapters), 
    color = "white") + scale_fill_viridis_c(option = "C", tran = "log", breaks = c(1, 
    3, 10, 40)) + theme_bw() + ggtitle("R-Lchapters per country")

Exercise 3: Raster and polygons

The tif file g100_06.tif contains land unse information for Europ. The file AX_Gebiet_Kreis.shp is a shape file of Baden WÃ¼rttemberg. Visualize the land use for BW. Do not forget to check crs. I have been trying, but I could not make the legend appear. Info for the legend in is the file clc_legend.csv. May be you find out!

# Load data

require(raster)
require(sf)

corine <- raster("data/g100_06.tif")  # loads the raster
shape <- st_read("data/AX_Gebiet_Kreis.shp")  # loads the vector data

## Reading layer `AX_Gebiet_Kreis' from data source `C:\Users\ichas\Documents\RLadies\geoR\rmarkdown\data\AX_Gebiet_Kreis.shp' using driver `ESRI Shapefile'
## Simple feature collection with 44 features and 2 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 388337.8 ymin: 5265159 xmax: 610069.4 ymax: 5515632
## epsg (SRID):    NA
## proj4string:    +proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs

# Inspect and plot the data Different functions to see the info that is
# encoded in our shape
dplyr::glimpse(shape)

## Observations: 44
## Variables: 3
## $ Name      <fct> LÃ¶rrach, Waldshut, Bodenseekreis, Ravensburg, Breisg...
## $ SchlÃ¼ssel <fct> 08336, 08337, 08435, 08436, 08315, 08326, 08327, 083...
## $ geometry  <MULTIPOLYGON [m]> MULTIPOLYGON (((418041.5 53..., MULTIPO...

colnames(shape)

## [1] "Name"      "SchlÃ¼ssel" "geometry"

summary(shape)

##               Name      SchlÃ¼ssel           geometry 
##  Heilbronn      : 2   08111  : 1   MULTIPOLYGON :44  
##  Karlsruhe      : 2   08115  : 1   epsg:NA      : 0  
##  Alb-Donau-Kreis: 1   08116  : 1   +proj=utm ...: 0  
##  Baden-Baden    : 1   08117  : 1                     
##  Biberach       : 1   08118  : 1                     
##  BÃ¶blingen      : 1   08119  : 1                     
##  (Other)        :36   (Other):38

plot(shape)

# Plot the raster

plot(corine)

# Check coordinate reference system The function to get the coor. system is
# different between shapes and rasters
st_crs(shape)  #get the projection

## Coordinate Reference System:
##   No EPSG code
##   proj4string: "+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs"

crs(corine)  #get the projection

## CRS arguments:
##  +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000
## +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

st_crs(shape)$proj4string == crs(corine, asText = T)

## [1] FALSE

shp.Transformed <- st_transform(shape, crs = crs(corine, asText = T))  #transform coordinate sistem

plot(st_geometry(shp.Transformed))  # plots only the geometry and not the associated values

st_crs(shp.Transformed)$proj4string == crs(corine, asText = T)  #check again if they are the same coordinate sistem

## [1] TRUE

# For extraction the shape file has to be of class 'Spatial Polygons Data
# Frame'
shp.Transformed.sp <- as(shp.Transformed, "Spatial")
class(shp.Transformed.sp)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

We get the extent (the coordinates of the extreme points) for our shape file of Baden WÃ¼rttemberg.

# get extent of Baden-WÃ¼rttemberg
corine.shp.extent <- extent(shp.Transformed.sp)
corine.shp.extent

## class       : Extent 
## xmin        : 4134156 
## xmax        : 4357500 
## ymin        : 2715907 
## ymax        : 2964364

# We crop the land use data by the extend of BW. Crop funciton from the
# raster package
corine.shp <- crop(corine, corine.shp.extent, snap = "out")
plot(corine.shp)

# mask BW - takes about 10 minutes
corine.shp <- mask(corine.shp, shp.Transformed.sp)

# 7) save corine.bav as a new file
writeRaster(corine.shp, filename = "data/corineBW.tif", format = "GTiff", overwrite = TRUE)
# loads it again
corine.shp <- raster("data/corineBW.tif")
plot(corine.shp, main = "Land Use")
scalebar(10, xy = click(), type = "line")

require(tmap)

# Visualize 1) import legend from a csv with readr package
legend <- readr::read_csv2("data/clc_legend.csv")
# note a special feature: attribute progress
`?`(readr::read_csv2)
# how does the legend look like?
knitr::kable(head(legend))

GRID_CODE	CLC_CODE	LABEL1	LABEL2	LABEL3	RGB
1	111	Artificial surfaces	Urban fabric	Continuous urban fabric	230-000-077
2	112	Artificial surfaces	Urban fabric	Discontinuous urban fabric	255-000-000
3	121	Artificial surfaces	Industrial, commercial and transport units	Industrial or commercial units	204-077-242
4	122	Artificial surfaces	Industrial, commercial and transport units	Road and rail networks and associated land	204-000-000
5	123	Artificial surfaces	Industrial, commercial and transport units	Port areas	230-204-204
6	124	Artificial surfaces	Industrial, commercial and transport units	Airports	230-204-230

## 4) plot map using tmap package

tmap_mode("plot")

tm_shape(corine.shp) + tm_raster(style = "cat", title = "Categories") + tm_layout(legend.position = c("left", 
    "bottom")) + tm_shape(shape) + tm_borders()

Exercise 4: Comparing bird observations in summer and winter

Here we will start justo from a data frame. The df has info of bird observations of one species (Pyrocephalus rubinus) in South and Central America. This bird is migratory, we want to visualize the location of the observations in different moments of the year. How would you do it?

The df is quite big, so to speed up I woud suggest to use recent observatins, e.g. after 2016.

require(lubridate)
require(ggplot2)
require(sf)
require(ggspatial)
require(dplyr)
require(tidyr)
require(rnaturalearth)
require(rnaturalearthdata)


load("data/Pyrocephalus-rubinus.RData")


Py_rubunis <- birds_df %>% mutate(m = month(observation_date), y = year(observation_date)) %>% 
    dplyr::select(m, y, latitude, longitude, species_observed) %>% dplyr::filter(!is.na(latitude)) %>% 
    dplyr::filter(!is.na(longitude)) %>% dplyr::filter(m %in% c(1, 7)) %>% dplyr::filter(species_observed == 
    TRUE) %>% dplyr::filter(y > 2016)



# We get
crs <- "+proj=longlat +datum=WGS84 +no_defs"  #EPSG 4326
Py_rubunis_geo <- st_as_sf(Py_rubunis, coords = c("longitude", "latitude"), 
    crs = crs)

world <- ne_countries(scale = "medium", returnclass = "sf")

Option 1 - Using ggplot and geom_sf

ggplot(data = world) + geom_sf() + theme_bw() + geom_sf(data = Py_rubunis_geo, 
    aes(color = as.factor(m))) + xlab("Longitude") + ylab("Latitude") + coord_sf(xlim = c(-120.15, 
    -25), ylim = c(40, -75)) + ggtitle("Observations of P. rubinus in January and July")

Option 2 - Using tmap

load("data/Pyrocephalus-rubinus.RData")


require(tmap)
require(sf)
require(rnaturalearth)
require(rnaturalearthdata)
require(dplyr)
require(tidyr)
require(lubridate)

world <- ne_countries(scale = "medium", returnclass = "sf")  # this is another function to get polygons of countries. 

countries <- c("Argentina", "Bolivia", "Brazil", "Chile", "Colombia", "Ecuador", 
    "Guyana", "Paraguay", "Peru", "Suriname", "Uruguay", "Venezuela", "Mexico", 
    "Belize", "Guatemala", "El Salvador", "Honduras", "Nicaragua", "Costa Rica", 
    "Panama")

sc_america <- world %>% dplyr::filter(name %in% countries)


Py_rubunis <- birds_df %>% mutate(m = month(observation_date), y = year(observation_date)) %>% 
    dplyr::select(m, y, latitude, longitude, species_observed) %>% dplyr::filter(!is.na(latitude)) %>% 
    dplyr::filter(!is.na(longitude)) %>% dplyr::filter(species_observed == TRUE) %>% 
    dplyr::filter(y > 2016)

crs <- "+proj=longlat +datum=WGS84 +no_defs"  #EPSG 4326

Py_rubunis_geo <- st_as_sf(Py_rubunis, coords = c("longitude", "latitude"), 
    crs = crs)


Py_rubunis_1 <- Py_rubunis %>% dplyr::filter(m %in% c(1))
Py_rubunis_4 <- Py_rubunis %>% dplyr::filter(m %in% c(4))
Py_rubunis_7 <- Py_rubunis %>% dplyr::filter(m %in% c(7))
Py_rubunis_10 <- Py_rubunis %>% dplyr::filter(m %in% c(10))

Py_rubunis_geo_1 <- st_as_sf(Py_rubunis_1, coords = c("longitude", "latitude"), 
    crs = crs)
Py_rubunis_geo_4 <- st_as_sf(Py_rubunis_4, coords = c("longitude", "latitude"), 
    crs = crs)
Py_rubunis_geo_7 <- st_as_sf(Py_rubunis_7, coords = c("longitude", "latitude"), 
    crs = crs)
Py_rubunis_geo_10 <- st_as_sf(Py_rubunis_10, coords = c("longitude", "latitude"), 
    crs = crs)


tmap_mode("view")

tm_shape(sc_america) + tm_borders() + tm_polygons() + tm_shape(Py_rubunis_geo_1) + 
    tm_symbols(col = "blue", scale = 0.5) + tm_shape(Py_rubunis_geo_4) + tm_symbols(col = "green", 
    scale = 0.5) + tm_shape(Py_rubunis_geo_7) + tm_symbols(col = "orange", scale = 0.5) + 
    tm_shape(Py_rubunis_geo_10) + tm_symbols(col = "red", scale = 0.5)

# tm_compass() + tm_scale_bar() + tm_legend(show = FALSE)

Visualizing Geo-Data in R

Elisa Schneider

5 July 2019

Difference between raster and vector representations

Simplest example

Information already availabe in R

From a couple of countries to the whole world: Plotting simple and nice maps of the world

1. We can just create a map from the world. We will do this using the library `ggplot2`.

2. We can change the colors, add text …

Coordinate reference systems (crs). Going from the rear world (the earth) to the simplifyed model (the map)

After choosing elipsoid and datum we have to project from 3D to 2D.

Make Choropleth Map

1. Load Data

2. Tidy the data set to have one column with country and another with life expectancy

3. Get data of the polygons of each country and merge map and life expectancy data

4. Plot with ggpoligon

5. Add points

6. Add scale and north

Rasters and polygons

1. Load required libraries

2. Load the raster files and do some calculations

This are some usefull links to find data:

3. Load the shape files

4. Check reference system

5. Crop the raster to the shape file extent

6. Plot with `tmap`

Hands on work!

Exercise 1: Geo-information available in R.

Option 1 - Using only ggplot2

Option 2 - Usinng sf objects

Option 3 tmap

Exercise 2: Geo-information from a CSV. Investigate the data frame. What can we learn from the data? which map could we do?

Option 1 – Only ggplot

Exercise 3: Raster and polygons

Exercise 4: Comparing bird observations in summer and winter

Option 1 - Using ggplot and geom_sf

Option 2 - Using tmap

Visualizing Geo-Data in R

Elisa Schneider

5 July 2019

Difference between raster and vector representations

Simplest example

Information already availabe in R

From a couple of countries to the whole world: Plotting simple and nice maps of the world

1. We can just create a map from the world. We will do this using the library ggplot2.

2. We can change the colors, add text …

Coordinate reference systems (crs). Going from the rear world (the earth) to the simplifyed model (the map)

After choosing elipsoid and datum we have to project from 3D to 2D.

Make Choropleth Map

1. Load Data

2. Tidy the data set to have one column with country and another with life expectancy

3. Get data of the polygons of each country and merge map and life expectancy data

4. Plot with ggpoligon

5. Add points

6. Add scale and north

Rasters and polygons

1. Load required libraries

2. Load the raster files and do some calculations

This are some usefull links to find data:

3. Load the shape files

4. Check reference system

5. Crop the raster to the shape file extent

6. Plot with tmap

Hands on work!

Exercise 1: Geo-information available in R.

Option 1 - Using only ggplot2

Option 2 - Usinng sf objects

Option 3 tmap

Exercise 2: Geo-information from a CSV. Investigate the data frame. What can we learn from the data? which map could we do?

Option 1 – Only ggplot

Exercise 3: Raster and polygons

Exercise 4: Comparing bird observations in summer and winter

Option 1 - Using ggplot and geom_sf

Option 2 - Using tmap

1. We can just create a map from the world. We will do this using the library `ggplot2`.

6. Plot with `tmap`